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Abstract 

In this paper we consider the numerical approximation of solutions to 
linear functional differential equations using the so-called Legendre-tau 
method. The functional differential equation is first reformulated as a 
partial differential equation with a non-local boundary condition involving 
time-differentiation. The approximate solution is then represented as a 
truncated Legendre series with time varying coefficients which satisfy a 
certain system of ordinary differential equations. The method is very easy to 
code and yields very accurate approximations. Convergence is established, 
various numerical examples are presented, and comparison between the latter 
and cubic spline approximations is made. 
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Introduction 


In this paper we consider Legendre-tau approximations of solutions of 
functional differential equations (FDEs). The tau method, invented by Lanczos 
in 1938 [12], is one of several approximation techniques which are referred to 
as a spectral method [9] . Spectral methods have been used in numerical 
computations for a wide class of partial differential equations (PDEs). In 
this paper, we view the original FDE as a PDE u fc = u x with boundary 
conditions which involve time differentiation. The Legendre-tau method is 
based upon representing the approximate solution as a truncated series of 
Legendre polynomials. The evolution equation for the expansion coefficients 
is then determined by substituting the series into the above PDE and by 
imposing the boundary conditions. To our knowledge, the use of the tau method 
in approximating solutions of FDEs and the specific manner in which it is used 
(i.e., applying the tau method to a reformulated PDE with boundary conditions 
involving time evolution) are new. The idea of formulating FDEs as Cauchy 
problems on an appropriate Hilbert space is not new. Within this framework. 
Banks and Kappel [2] make use of approximation results from linear semigroup 
theory (in particular, the Trotter-Kato theorem) to establish the convergence 
of numerical schemes based upon splines. 

In this paper our considerations are restricted to linear autonomous FDEs 
of retarded type. The use of our ideas for (i) nonautonomous , (ii) nonlinear, 
(iii) neutral-type or (iv) integro-dif ferential systems will be discussed 
elsewhere. We add here that the tau method should prove useful in dealing 
with partial differential equations with boundary conditions which involve 
time differentiation. Our main goal is the application of the tau method to 
optimal control and parameter estimation problems. As will be discussed in 
section 7 the tau method may offer considerable improvements over other 
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methods (e.g.- those discussed in [2], [3]) in many instances. One reason for 
this is that the semigroup (S(t) : t > 0} associated with a retarded FDE has 
the property that the range of S(t) is contained in V(A^) for each t > kr 
where A is the infinitesimal generator of (S(t) : t > 0} with domain V(A) 
and r is the longest delay time appearing in the FDE. Thus, the regularity 
of solutions increases with time. In section 3, we will see that, in such a 
case, approximations by orthogonal polynomials are quite powerful. 

The following is a brief summary of the contents of this paper. In 
section 2, we review the equivalence results between FDEs and abstract Cauchy 
problems on the product space if 1 * Lj. In section 3, we recall various 
properties of Legendre polynomials including certain estimates which are 
needed to establish convergence. Section 4 is concerned with the developemnt 
of the numerical scheme based upon the Legendre-tau approximation for solving 
the class of FDEs under consideration. In order to establish the numerical 
convergence of such approximations to the actual solution in the IFf 1 x 
norm, we first show, in section 5, that convergence holds in a stronger norm 
under the special assumptions that the initial data is sufficiently regular 

and that inhomogeneous forcing term f is identically zero. In section 6, we 

lo c 

then extend our result to the inhomogeneous case wherein f e and show 

that the sequence obtained by truncating the last term in each of the 
Legendre-tau approximations converges to the actual solution in the S' x 
norm whenever the initial data lies in if 1 x L 2 » Finally, in section 7, we 
present numerical results and compare these results for cubic spline 
approximations discussed in [2]. 

Throughout this paper the following notation will be used, r > 0 stands 
for the longest delay time appearing in the FDE. The Hilbert space of IP- 
valued square integrable functions on the interval [a.b] is denoted by 
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[a,b] ; if 1 ) • When the underlying space and interval can be understood from 
the context, we will abbreviate the notation and simply write l^* 

L2° C ([0,°°); I?), or L^ 00 , is the space of IP-valued locally square 
integrable functions on the semi-infinite interval [0,°°). is the 

Sobolev space of IP-valued functions f on a compact interval with f(k-l) 

(k) 

absolutely continuous and f v e L^. The Banach space of IP-valued 
continuous functions on the interval [-r,0] is denoted by C. We denote 
by Z the product space IP * ( [-r,0] ; IP) . Given an element 

z e Z, n e IP and $ e denote the two coordinates of z : z *• (n <}>).The 
bracket <*,*> stands for the inner product in the Hilbert space H, and the 
subscript for the underlying Hilbert space will be omitted when understood 
from the context. II • II denotes the norm for elements of a Banach space and 
for operators between Banach spaces, while | • | denotes the Euclidean norm in 
IP. 

If X and Y are Banach spaces, then the space of bounded operators 
from X to Y is denoted by L(X,Y). V (A) denotes the domain of a linear 
operator A • Xj denotes the characteristic function of the interval I. 

Given a measurable function x : [-r,°°) ♦TP and t > 0, the function 
x fc : [-r,0] ->-xP is defined by x t (9) = x(t+9), 9 e [-r,0] . Finally, for 
any function <}> of the independent variable 9, we shall use either <j> or 
to denote the derivative of with respect to 9. 

2. Linear Hereditary Differential Equations 

In this section, we state the type of equations to be considered and 
recall some results for such equations which are important for the discussion 


to follow. 
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Given 

equation 


(n,<}>) e Z and f e [0,°°); IP) , we consider the functional 


t t 

x(t) = n + / Dx ds + / f(s)ds 


0 


0 


x 0 = ** 


where D : V(D) C ( [-r , 0] ; ]RP) iP has the form 


( 2 . 1 ) 


,0 

D <J> ■= / dp(6)<K 0 ) 

-r 


( 2 . 2 ) 


with P a matrix-valued function of bounded variation on [-r,0] . As an 
example, consider 

P(9) - ). A, x, 0 rtl (0)+j A(s) ds . (2.3) 


A X/ 0 ni ( 0 ) + / A(s) ds . 
i=0 1 { i’ UJ 


0 

* 

-r 


where -r = 0 <•••< 0_ = 0 and A* and A(*) are n x n matrices, the 

m 0 t 

elements of the latter being integrable on [-r,0]. Then 


Dx = ? A x(t+0.) + J° A(0)x(t+0)d9, 
i=0 1 1 -r 


t > 0. 


(2.4) 


1 OP 

For (n , <}>) e Z and f £ , there exists a unique solution x to 


(2.1) in L 2 ( I-r,T] ; K 1 ) for any T > 0 such that 


z(t) = (x(t),x ) e Z, 


(2.5) 


is a Z-valued continuous function which depends continuously on (n,<)>) e Z 
and f e L, 2 ° C for each t > 0 (see [3], [5], [8], e.g.). Note that x(*) 
is absolutely continuous on [0,T] and satisfies the differential version of 
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equation (2.1): 

(t) = Dx t + f(t), t > 0 (2.6) 



For t > 0, define S(t) : Z ■* Z by S (t ) (t) » 4> ) ■ (x(t),x fc ) where x is 
the homogeneous solution of (2.1) (i.e., f = 0). Then (S(t) : t > 0} forms 
a strongly continuous semigroup on Z. The following results are now standard 
if one deals with FDEs in the state space Z ([3], [4], [15]). 

Lemma 2.1: 

(i) ^ A denotes the infinitesimal generator of (S(t) : t > 0}, then 

V(A) » { (n , 4> ) e Z : } e l 2 and n ■ <J>(0)} 
and for (<HO),<f0 e p(A), 

A ( 4> (0) ,<p) = (D<M). 

(ii) If z(0) = (<K0),<J>) E V(A), and f e L^ oc , then 

z(t) = S(t)z(0) + / S(t-s) f(s)ds satisfies the equation 

0 

(t) - Az(t) + f(t), t > 0, (2.7) 

in Z where 8 : If 1 Z is defined by 

8x ■ (x,0) e Z for x e IP. 
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(iii) For k > 2, 

V (A k ) C { (4>(0) ,<{)) £ Z : $(0) = D* and <|> £ H k } , 
and is dense in Z. 


3. Properties of Legendre Polynomials 

In this section, we review some properties of Legendre polynomials ([11], 
[13], e.g.). 

The Legendre polynomial of degree n, p n (x), -1 < x < 1, can be defined 
as the solution of the differential equation 


d^ ( (1 " x2) % (x) ) + "(n+l)p(x) - 0, 


(3.1) 


which satisfies p ( 1 ) = 1. Thus, p Q (x) = 1, pj(x) « x, p 2 (x) =1 ^(3x 2 -l), and 
so on. The Legendre polynomials ^P^kX) sat ^ s ^y the orthogonality relation 


1 

/ p (x)p (x)dx 
, m n 


_2 

2n+l 



(3.2) 


and they form a basis for L 2 (-l,l) : any f e L 2 (-l,l) can be written as 


f = £ f p 

n>0 n n 


where 


2n+l 


"n 2 


/ f(x)p (x)dx 


-1 


n 


l,f,l L " l TTT f2 * 
L 2 n>0 2n+1 n 


with 
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They possess the recursion formula 

(n+l)p . (x) = (2n+l)xp (x) - np .(x). (3.3) 

n+i n n-i 


From this, we have 


P n (±l) = (±l) n , lp n (x)l < 1, Ix| < 1, 

and 

p (±1) = (±l) n n(n+l)/2. 
n 


If f is represented as 


then 


N 

I 

n-=0 


f n P_ > 
n n 


f *» l b p where b = (2n+l ) l f, • 
a n n n * . i 

n=0 k**n+l 

k+n odd 


(3.4) 


For any positive integer N, let P N be the orthogonal projection of 

i tN 

l >2 onto the subspace spanned by '■Pk^k=0* Then we have the following error 
estimates [6]: 

Lemma 3.1: For any real s > 0 , there exists a constant K such that 


Ilf - P N f ll T < KN _S II f II 

L 2 H 8 


Lemma 3.2: For any real s and a such that 1 < s < 0, there exists 


a constant K such that 
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Ilf - P N f II < KN 2s “ a “^ 2 llflt 

H S H° 


The next lemma gives an error estimate in the sup remum norm: 


Lemma 3.3: For any positive Integer m there exists a constant K 

such that 

|f(x) - P N f(x)| < KN“ 2mfl IIfll , . 

H 2m 

Proof : Let us denote by L the differential operator 

(L£)(X) 5 W)' 

From (3.1), we have 


1 l 1 

/ f(x)p n (x)d = ( n T n + i ~ )^ m ^ ( L m p n ) (x)f (x)dx, 


n > 1. 


Since L is symmetric, 


1 1 1 
/ f(x)p n (x)dx = ( n - (n+ ' i) ) m / ( L m f) (x)p n (x)dx, 


n 


> 1. 


It follows that if (f } are the Legendre coefficients of f, then 

n n>u 


£ „ ■ ferry)"’ [ t < i"£><i>p”<*>dx 5 ( 


1 im 


n(n+l) j g n’ 


where {g n ) n> Q are coefficients of L m f- Thus, for any M > N 
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and |x| < 1 , 


MM M 

IP f(x) - Pf(x)| - | l fp (x) I 


n=N+l 
M 

l 

n=N+l 


n n 


< l If I = l lg n l 

n=N+l n k=N+l n(n+1) n 
f ( y ( 1 «.2m 2n+l-t V? r r 2 , ,2^ 

< 1 1 tery 5 —> *• _ - . , 2n+r isj 3 


n=N+l 


n=N+l 




N 


f». 


(4m4.2)- l *(-i f rT)"- 1/ 2| L-*. . 


Here we have used Ip^Cx) I < 1 for |x| <1 and 


r 1 i2m 2n+l < r ( 1 2m 2x+l 

"n(n+l ) ' 2 * } . l x(x+l) J 2 

n-i 


dx. 


It now follows that {P N f)„^„ Is a Cauchy sequence In C[-l,l] and hence 

N>u 

that P N f converges uniformly to f. Letting M “ above, we obtain 


|f(x) - P N f(x)| < (4m + 2)" 1/2 (jiTjj+r)” i mfll L • 


Since L m is a differential operator of order 2m with continuous 

coefficients on [-1,1], there exists a constant C m for each m > 0 such 
that 

n L m f« < c iif ii 0 , 

m 


which completes the proof 
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oo N * 

If f is C , then, from Lemma 3.1, the error Ilf-P fH decreases more 

rapidly than any power of 1/N. This is usually referred to as "infinite 

order" approximation. 


4. Legendre-Tau Approximation 

In this section, we discuss the Legendre-tau approximation of solutions 
to (2.1). Throughout the rest of the paper, it is convenient to modify the 
Legendre polynomials by shifting their fundamental interval from [—1,1] to 
[-2,0]. Thus, the p^'s will be understood to be orthogonal on [-2,0] with 
P k (-2) - (-l) k and p k (0) = 1. 

For simplicity of exposition, we assume that r = 2. See the remark at 
the end of the section for the general case. If z(t.9) = x(t+9), 9 e [-2,0], 
where x is the solution to (2.1) with initial data (n,<{>) e I XA). then, 
according to Lemma 2.1, z satisfies 

If (t ’ 0) = 41 (t ’ 9) ’ 9 e [-2,0] (4.1) 

dz 

(t,0) = / du(0)j5(t,9) + f(t). (4.2) 

dt -2 

N 

The approximate solution z (t,9) is assumed to be expanded in a Legendre 
series: 


N 

z 


(t.9) 


I a k (t)p k (9). 


a k e IP. (4.3) 


The tau approximation [9] of (4.1) and (4.2) is as follows. Note that from 
(3.3), 
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N N-l 


where 

b (t) = (2k+l) 5> a (t) . (4.4) 

j=k+l J 
j+k odd 


3 „N 

Equating 


with 



leads to the N equations 


d_ 

dt 


a k (t) 


b k (t) 


0 < k < N-l. (4.5) 


The essence of the tau method Is that the boundary condition (4.2) is then 
Imposed to determine an equation for a N - From (4.2) we obtain 

. N 0 

7 - ( l a, (t) ) “ / du( 6 )z (t,9) + f(t), 
k =0 -2 

or 

, N-l 0 „ 

7 - a„( t) = - l b, (t) + / dh(0)z N (t+0) + f(t). (4.6) 

k =0 -2 

Hence, from (4.3) - (4.5), we obtain a system of ordinary differential 
equations for aQ,***,a^: 


d_ 

dt 


a N (t) 


A N a N (t) + B N f(t) , 


(4.7) 


where 


col(a Q 


,a N ) and 


B 


N 


© 


I, 


where 


e^ e is given by ^ = col( 0 , 


0 , 1 ), I is the n*n identity 
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matrix, and © denotes Kronecker product. For the case where N is even; 


is given by 


where 


.N _ ,N N 

A = A-, + A , 

0 p 


.N 


0 

0 

0 


1 

0 

0 


0 

-1 


0 

3 

0 


0 

-3 


1 

0 

5 


0 

-6 


0 

3 

0 


0 

-10 


1 0 

0 3 

5 0 


2N-3 

0 


0 

2N-1 

N (N+l ) 
2 . 


© I (4.8) 


and 


a n = 

P 


,N 




(4.9) 


with 


w 0 

D k = ^ d ^ 9 >P k > 


0 < k < N. 


(4.10) 


Let us introduce the projection operator L 1 * on Z. For any 
z - (n,<J>) e Z, L N is defined by 


l N z = (n, f a.pj, 
k=0 K K 


(4.11) 
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where, for 0 < k < N-l, 


2k+l 
l k 2 


f <j>(0)p. (8)d6, 


-2 


and 


N-l N-l 

a N = n " ^ a k P k (0) = n “ ^ V 

k=0 K k=0 


The tau method can then be interpreted as follows. Let 


N , N N s 

z (t) «(z"(t,0), z W (t,*)) e Z 


where 


N 


Z N (t.0) = l a k (t)p k (0) 
k=0 


as in (4.3). Then the approximate solution z N (t) satisfies 


d N,_ T N . N . . . t n c 

z (t) = L A z (t) + L 8f(t), 


(4.12) 


in Z. Indeed, 


with 


N 

L N Az N (t) + L N Bf(t) = (Dz N (t,*) + f(t), I b. (t)p. ) e Z, 

k=0 k 


m N-l 

b N (t) = Dz W (t,-) + f(t) - l b k (t) 

k=0 


= D z N (t , • ) + f(t) - z N (t,0). 


Thus, the first component of (4.12) is equivalent to (4.6) and the second 
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component is equivalent to (4.5). Alternatively, the approximate solution 
z N (t) is the exact solution to the modified equation 


z N (t) = A z N (t) + Bf(t) +t N (t). (4.13) 

in Z where T N (t) = (o,b N (t)p N ) e Z. 


Remark : For the general case, l.e., r * 2, the matrix 

(4.7) is replaced by 


a" - i A? + A N 
r 0 u 


aN 


where 


-N 


0 

I 


-r 


dU (0 ) P j c (~~) , 


0 < k < N, 


appearing in 


N 

replaces in (4.9). In particular, if the delay operator D is given by 

(2.3) and (2.4), then has the form 


-N 


m 

I 

i=0 


29. 


A i P k 


(“A 


0 

+ / 

-r 


A(0) Pk^ d0 » 


for 0 < k < N. 


5. Convergence Proof 

In this section, we show that the approximate solution z N converges to 
the true solution z of (2.1) in the case where the forcing term of f = 0 
and the initial data lies in 0(A). The cases where f * 0 and the initial 
data lies merely in Z will be dealt with in the following section. The 
principal result is based upon the Trotter-Kato theorem (see Theorem 4.6 in 
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[14]). 

Theorem 5.1: Let S(t) and S N (t), N > 1, Jse^ Cn -semigroups acting on a 
Banach apace X with Infinitesimal generators A and respectively. 
Assume that the following conditions are satisfied: 


(1) (stability) There exists a constant w such that A - “I and 
N 

A “ ul are dissipative on X. 

(11) (consistency) There exists a subset V contained In V(A ) 

0 ° 

D 0 v { A ) which together with (Al - A)V for some A > 0 Is dense In X 
N -! 

and such that A <f> + A^ for all 4> e p. 

Then for all <)> e X, II S N (t)<J>-S (t ) 4> n -*■ 0 uniformly on bounded t-intervals . 

We begin by introducing the Hilbert space X = H^[~r,0] with norm 

2 2 r® • 2 

*♦*1 s |«KO)r + I M>(0)| d0 

-r 


and Inner-product 


r° • 

<4>,ip> 1 = «J>(0),iK0)> I? i + / <H s ),^( 9 )>ifd9. 

-r 

It is readily established that X is equivalent to V(A) C Z (see Lemma 2.1 

0 0 0 

for definition of V(A) ) with graph norm ( ,lzl, p(A) E ' lzl1 + "Ad there 

exist constants 0 < c < C < 00 such that 


d + lj < «(*(0) f *l ( } < C I <f> II x 


(5.1) 
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Let us denote this isomorphism between X and V(A) by E : V (A ) ■* X : 

E(M0),<j>) = M Since A generates a CQ-semigroup on V(A) C Z with graph 
norm, EAE - ^ generates a CQ-semigroup, which we denote by { S ( t ) : t>0}, on 
X. Moreover, it is easily seen that 

EAE = 30 

and 

p(EAE _1 ) = !♦ ex : i e X and MO) = D<|> } 


For the sake of convenience, however, we will use the same symbols to denote 
corresponding elements in and operators on V(A) and X. except in one 
instance where greater clarity is required. Thus, on X. we have 


A = 


3 _ 

30 


and 


PCA) = {<f> e X : 1 e X and MO) = D$} 


Without loss of generality, we can assume that r = 2. For any positive 
integer N, we define the finite dimensional subspace X^C X by 


N 


X N = {<}> e X : M0) - 1 a, p, (6 ) , a,e If 1 } 

k=0 K K K 


and define the orthogonal projection Q N of X onto X N by 


CQ N 4>) (0) = MO) 


(5.2) 


d_N , _ pN-1 1 
<P = P 9 


(5.3) 


In fact, if 


q n * 


- ? 


k=0 


b k P k* 


then 
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d-N. 

1ft 41 = 


N-l 

l b k P k * 
k=0 


r l N 

so that, from (3.4), are uniquely determined by (5.2). bQ Is then 

determined by (5.1): 


Cq N 4>) (0) - l b. = <f>(0) 

k=0 K 


N 

*(0) - I b. 
k=l K 


Let us define A X -+■ by 


N N -1 N 

A <J> = EL AE Q <j> for $ e X 


(5.4) 


where A Is understood as In Lemma 2.1 and L N is defined by (4.11). 
Returning to our above-mentioned convention, (4.12) can then be written as 

^z N (t) = A N z N (t) + B N f(t) in X (5.5) 

where, for N > 1 , B N = L N B» 

Although it is our goal in this section to merely prove convergence for 
the case where f = 0, it is simiple enough to include, in the latter part of 
the following Lemma, a statement concerning the more general inhomogeneous 
case which will be used in the next section. 


Lemma 5 . 2 (stability) There exists a constant m e !R such that, for all 
N 

N > 1, A - ml is dissipative on X, i .e . 


.N , , 

<A 4> > < 


con^n 


for 4> e x 
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Moreover, if z^(t) is the solution to (5.5) with z“(0) = for <{> e X, then 




lz N (t) II J < e 2ut llz N (0)ll 2 + / t e 2aj(t_s) |f(s) | 2 ds (5.6) 


Proof : For any e X^ and f e I?, it follows from (5.4) and (4.13) that 


0 


^ N A N f,<f> N > 1 = <D<j> N +f,Ao)> + / 4 N (0)+b N p N (6)^ N (0)>d0 (5.7) 


where 


-2 


N *N 

b N = d<j> n + f - Ao) 


(5.8) 


Let I denote the lntergral term in (5.7). Then 


1 = ^2 fe^ N ( 0 )l 2de + / 2 <b N ,?( 0 )>P N ( 0 )d0 


1 *N 2 1 . *N 2 *N N # N 

ji-j) (o)r - ji<t> (-2)r + <b N ,<j> (0)-(-)V(-2)> 


,0 "N 

+ / 2 <r(9),b N >p N (9)d0 


(5.9) 


where we have used pjj(O) = 1 and pjj(-2) = (-1) N . Note that the last term 

'*N __ 

on the right-hand side vanishes because <j> is a vector in R 1 whose 
elements are polynomials of degree less than N and hence orthogonal to 
Pjj. Thus, from (5.9), 


1 •tj 9 1 *N .2 N *N *N N*N 

I --§|4> (0)| Z - j\V (-2)1 + <Df+H (0),H0)-(-)V(-2)> 

1 *N 2 N *N *N l.*N .2 

--■jir(o)r+ (-) N «f> N (0) , <t> N (-2)> (~2)r 
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N *N N»N 1 N 2 

+ <D<j> +f,4> (O)-(-) * (-2)> < 2 ID* +f| 

2 2 

where we have used the inequality 2 <x,y>j£i < |x| + |y| • From this and 

(5.7), we find 

NNNN NN 1N2 

< a % +b f,* >j < <D<Af,Ao)> +-jiD*“+fr 

NN NN 1 2 1 N 2 

= <D*“,r(o>> + <f ,D<p+<j> (0)> +i|fr + jiD*"r 

N N 1 , 2 N N 

- <D* ,* (0)> + { | f | +2<f ,Di{> +9 (0 )> 

+ |dA* N ( 0)| 2 } - ^ID4> N -H> N C0) | 2 -|lD4> N 1 2 

- jlfWVfO)! 2 ‘ lA°)l 2 < |d* n +* n ( 0)| 2 -+ |f I 2 

Note that D e /.(X, IFf 1 ) with 

|D<j> | < 01*1 (5.10) 

for some positive constant 3 < ". Thus 

<A N * N +B N f ,<)) N > 1 < toll <J> N ll^ + |f | 2 (5.11) 

with a) = (1+3) 2 

Since Q^L^z = L^z for all z e Z and is symmetric with respect to 


w e have, for all * e X, 
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N N N N N N N 

<A <M>J = <Q A Q <M > 1 = <A Q <J>,Q 


It then follows from (5.10) that 


<A N 4>,4>> 1 < < toll< * )ll i 


,N 


which implies that A - <oI is dissipative on X. 
From (5.5) and (5.11), 


1 d , N, . „ 2 .N N, ..X. . N. , . „ N, x ,,2 , ... .,2 

2 dt z ^ l = z (t)+B f(t),z (t)>j < (ol. ( t ) II ^ + | f ( t ) | 


A standard agrument using Gronwall's inequality then yields (5.6). 


(Q.E.D.) 


Lemma 5.3 (consistency) Let 


p k = {<f> e X: cf> e H k and 4>(0)=D<f>} 
k k 

Then for k > 2,p and (Xl-/\)p for X e H sufficiently large are dense 
in X and, for k > 5,A^4> A 1 )* .in. ^ for <(> e p k . 

Proof : Since 0 k densely contained in 0(A) for k > 2, the denseness of 

0 k follows from that of 0(A) • Since 0^ = 0(A) and, for XeH sufficiently 
large 


(XI- A) 0(A) = X, 
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( Al— A )2/ is dense in X for k > 2.- 

For 4> e P(A) and Q^<}> = 4>^, 


(A N -AH = - £ + (d^. n 4 n (0))p n 


(5. 


in X for <p e IXA)* From (5.10), the definition of Q 1 * and Lemmas 3. 
3.3, 


*N * _9 rn 4- 1 • 

|4> C0)-4>C0) I < KN n<i>i 


,2m 


, ?N • , T _ -2mf 3 /2 „ • 

!l<j) -<}>ll < KN ll(}) II 


H 


2m 


and 


|d(^) N - 4| < (l+3)KN _2m ll 4>ll 


2m 


Since $(0) «= D<() for <}> e p(A), 


N *N -2m+l • 

|D<}> -<J> (0) I < (1+3)KN n (}>II 


H 


2m 


Using the formula (3.4), it is easy to verify that II = /N(N+1) . It 
follows from (5.12) that 


n Ca N -a)^u l < 4 N -4>n 1 + |d<} n -4) N (0) I(i+/N(N+1)) < (2+3)KN _2mf2 4n 


H 


2m 


Hence, if <p e V K , k > 5, 


n (A N -A)<}>i! , < (2-t-0)KN -2 II 4>n c o 

H 


12 ) 
1 - 


now 


(Q.E.D.) 
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Remark : Using the same "completing the square" argument as the 

in the proof of Lemma 5.2, it is easily verified that A - wl 
on X as well. 

Combining Theorem 5.1 with Lemmas 5.2 and 5.3, we have 

f N l 

Theorem 5.4: {S (t):t>0| denotes the semigroup on X 

A 1 *, N > 1, then, for all $ £ X, 

ns N (t)4>-s(t)<|)ii 1 -*• o 


one employed 
is disipative 


generated by 


uniformly on bounded t-intervals. 
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6. Convergence Proof (continued) 

In this section, we prove the convergence of our scheme for the cases 
where f ^ 0 and where the initial data (n,<f>) is merely assumed to lie in 
Z. Throughout this section, the operator A is defined as in Lemma 2.1. 

Let us first consider the case where f £ L^^tO, 00 ); if 1 ) and 
= 0. The solution to (2.1) is given by 

z(t) “ / S (t-s)Bf (s)ds (6.1) 

0 

Recall that for any f e if 1 , Bf ** (f,0) e Z. It is easily seen that 

A _1 Bf - (A -1 f,p 0 A -1 f) e z (6.2) 

0 

provided that A = / dy(0) is invertible. 

-r 

It is easily seen that the spectrum of A conicides with the zero set of 
the entire function 

0 xe 

c(X) = det[xi-/ e dp(6)] 

-r 

If c(0) = 0, i.e., if A is not invertible, we can choose any X e P(A), i.e., 
c(X) * 0, and consider y(t) = e ^ t x(t), x(t) being the solution to the 

initial-value problem (2.6). Then y(t) will he the solution to 

^jry(t) - -ly(t) + / e* 9 dp(0)y(t+0) + e * fc f(t) = D^y t + 6 
-r 

with initial data y(9) = e ^%(9), -r<6<0. 
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Clearly, this problem may be formulated as before on the product space 
Z; the corresponding generator A x will by construction satisfy 0 e p(A^): 

c A (w) = det[ ( co-t-X ) I — / e (<i)+X)0 dy(6)] 

-r 

c x (0) = det(A x ) = c(A) * 0 

Thus, without loss of generality, we may consider (6.2) above. 

Let us rewrite (6.1) as 

t 

z(t) = J AS (t-s)Cf (s)ds 
0 

where Cf = A -1 Bf. If f is continuously differentiable, It follows from 
[10, p.488] that 

z(t) = S(t)Cf(0) - Cf(t) + / S(t-s)Cf(s)ds (6.3) 

0 

From the definition of A N and B N » it is easy to show that 

(A N ) _1 8 N f = Cf for f e I?. 

(See Lemma 6.2 below). Applying the same argument to 

z N (t) = / A N S N (t-s)(A^) *Bf(s)ds 
0 

yields 

z N (t) = S N (t)Cf(0) - Cf(t) + / t S N (t-s)Cf(s)ds 

0 


(6.4) 
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Since Cf e P(A) for f e l? 1 , it then follows from Theoren 5.4 that 

S N (t)Cf S(t)Cf in X for f e TP 


and since dim(]EP) < °°, this implies 


ll(s N (t)-S(t))c« L(lPfX ) - 0 


Hence, from (6.3) and (6.4), if f is continuously differentiable. 


llz N (t)-z(t)H 1 + 0 


(6.5) 


uniformly on bounded t-intervals. 

According to the stability result of Lemma 5.2 (c.f. (5.6)), 



< 



( [0,t] ; ®P) 


( 6 . 6 ) 


According to the remark made prior to the statement of Theorem 5.4, the above 
estimate holds true with S N replaced by S: 


t 

1 / S(t-s)Bf (s)dsll 
0 


«t, 


fit 


L 2 ( [0, t] ; TP) . 


(6.7) 


Since the space of continuously differentiable functions on [0,t] is dense 
in L 2 (0,t) for all t > 0, a simple limit argument using (6.6) and (6.7) 
shows that (6.5) holds for any f e L^ 00 as well. 

Thus, combining the above with Theorem 5.4, we have 
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Theorem 6.1: For any Initial data (<l>(0),<|>) £ z with <p £ and 

f £ the approximate solution z^(t) to (4.12) converges strongly to 

z(t) _in X, uniformly on bounded t-lntervals . 

Note that Theorem 6.1 Is not applicable when the Initial data (n,<f>) Is 

merely assumed to lie in Z. Thus, to complete our study of convergence, we 

need to consider the case where ( ti , ) £ Z and f = 0. Let us define the 

N 

orthogonal projection II on Z by 

n N (n,<f ) - (n,P N-1 <|>) for (n,$) e Z. (6.8) 

Lemma 6.2: For z = (n,<f>) e Z, 

fl N r l„N _ -l n N 

(a J el z = EA Hz 

Proof : Recall that the isomorphism E:V(A) -*■ X and that 

A N = EL N AE _ ^Q N : X X N = Q N X; thus, there will in general be many solutions 
w £ V(A) to 

A N Ew - EL N z (6.9) 


We shall interpret (A N )”*EL N z to be the unique solution of (6.9) lying in 
X^. We must then show that a unique solution exists and has the indicated 
form. 


For z £ Z, let 


N— 1 


n N z = (n, l b fc p k ) = (n,<j>) 


k=0 


then 


N N-l 

L z = (n» l b, p, ) with bj. = n - I b = (n,<f>+(n-4(0))p ) 
k=0 K K k=0 k W 


n 
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N N N 

and it follows that L H z = L z. A short calculation shows that 
A -1 (h, = CKO),*) where 


-1 0 0 
iJ/(0) = A (n+D/ <j>(s)ds) - / <Ks)ds. 

6 


and thus, since <J> is a polynomial of degree N-l , * has degree N 
shows that 

_ ,-l _N „N 
E A n z e X 


N N 

Since Q :X -*■ X is a projection, 


aN^-IttN I7T N Ar.- 1 rt N r.A- 1 TT N 
A EA II z = EL AE 0 EA Hz 


m N n N 
EL n z 


= EL z, 


i.e., Ew = EA _1 n N z e X N is a solution of (6.9). 
To establish uniqueness, we must show that 


A N Ev ** 0 and Ev e => v =* 0 


Let Ev = 


N 

l Vk 

k=0 K K 


= <J> so that 


_-l . 
v = E <j> 


(*(0),*) e P (A) 


This 


Thus, A v = (D , 4* ) , and so 
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A N Ev = EL N AE -1 Q N Ev = EL N Av = e(d<}),^+(D<!>-^(0))p n ) = I + (D<f)-|(0))p N 

Since 4> is a polynomial of degree N-l, A^Ev = 0 and the orthogonality of 
the Legendre polynomials yields the fact that $ = <f>g = constant and 
D<J> 0 = 0. But 

D 4>o - D^» 0 p 0 = 

whence (c.f. (6.2)) <J>g = 0, or v = 0. 

(Q.E.D.) 

Theorem 6.3: For any z = (n,<(>) £ Z, 

N/ -1 N N 'v 

II n (e S (t)EL z - S ( t ) z ) ii z o 

uniformly on bounded t-intervals. 

Proof : Let 

d = d(t) = n N [E _ 1 S N (t)EL N z-S(t)z) 

thus 

d = n E As (t)(_A J EL z - n S(t)z 


It follows from Lemma 6.2 that 


d 


_N_-1. N_N . -1..N J n , . n Nf -li N . -hJ, -UN 

it e A s (t)E n z - n s(t)z = n [e a -Ae Js (t)E n z 

+ n N AE _1 (s N (t)-s(t))EA" 1 z + n N AE* 1 s N (t)EA" 1 (n N z-z) = dj-H +d 3 . 
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N 

Since S 


because 

lid 


by (5.1) 
intervals 

lid 


by Lemma 


and thus 
completes 


( t)E -1 n N z e X N , 

d l = n N (E -1 EL N AE _1 Q N -AE -1 )s N (t)EA" 1 n N z 
= n N (L N -i)AE _1 s N (t)EA" 1 n N z = o 

II^L^y = Il^y for any y £ Z. Since is orthogonal on Z, 

2 H Z < HAE“ 1 (s N (t)-S(t))EA" 1 zn z < IIE" 1 (s N (t)-S(t))EA~ 1 z» p(A) 

< Cll(s N (t)-S(t))EA _1 ztl l 

. Thus, by Theorem 5.4, Hd z (t) 0 uniformly on bounded t- 

. Likewise, 

„ 2 , 2 N , . 1 /- ,-N > 2 2 „_N. ...2 

3 n z < c us (t)EA (n z-zjiij^ < c • iis (t)“ L (x>* EA z-zJHj 

„ _2 2mt -2 -1 /--N \.2 

< C e *c HA (H 2-2 ^ 0(A) 

5.2 and the first inequality in (5.1). Thus 

I' d 3 H z < (f) 2e2 “ t { l| A _1 ( nN z-z) ll 2 +linNz - z!l z^ 

< (^-) 2 e 2wt * (uA _1 n^ (z) +i)nn N z-zii z 


lld 3 (t) H z 0 uniformly on bounced t-intervals as well. This 
the proof of Theorem 6.3. 
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7 . Numerical Results and Conclusions 

In this section we discuss some numerical examples which demonstrate the 
feasibility of the Legendre-tau approximation. For the purpose of comparison, 
we have also computed approximate solutions by using the cubic spline 
approximation (Sg) which is discussed in [2]. All computations were 
performed on a Control Data Corporation Cyber 170 model 730 at NASA Langley 
Research Center (LaRC) using software written in Fortran. The integration of 
the system of ordinary differential equations (ODEs) (4.6) was carried out by 
an IMSL routine (DVERK) employing the Runge-Kutta-Verner fifth and sixth order 
method. 

For the Legendre-tau approximation, implementation of the algorithms is 
almost as easy for the averaging approximation (AV) which is discussed in [3] 
and the first-order spline approximation (S^). In the following, we give the 
number of operations for each approximation to compute the right-hand side of 
(4.6) for the scalar system: 

4r x(t) = ax(t) + bx(t-r) + f(t). (7.1) 


Table I 


Method Number of Additions Number of Multiplications 


SP 

2N + 4 

N 

+ 

2 

AV 

N + 2 

N 

+ 

2 

S 1 

3N + 4 

3N 

+ 

5 

s 3 

11N + C x 

9N 

+ 

c 2 
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SP stands for the Legendre-tau approximation. For S 3 , Cj and C 2 are some 
positive constants (independent of N). Note that, for and S 3 , the 

operations which are required in order to perform the Cholesky decomposition 

N 

of (see, for the definition [2], p. 511) are not included. However these 
numbers may not reflect directly the CPU time required to solve the 
approximating system of ODEs. According to our calculations, for the same 
value of N, the Legendre-tau approximation is about four times as fast as 
S 3 . In addition, no storage space is necessary for the matrix A N appearing 
in (4.6) for system (7.1) because of its simple structure. For the general 
system (2.6) in if 1 , the last n rows of A n need to be stored in order to 
integrate the approximating system of ODEs using the DVERK routine. 

As will be evident from the numerical results for the examples presented 
below, both approximation methods (SP and S 3 ) behave about the same 
initially. The typical feature of the Legendre-tau method is that the 
relative error decreases with increasing time. This is not unexpected due to 
two facts: (i) the regularity of the solution to a problem which has smooth 
inhoragoeneous terms increases in time (as pointed out in the Introduction) and 
(ii) Lemma 3.1. For all examples, on the interval where the solution is 
infinitely differentiable, the rate of convergence seems to be infinite order. 
In contrast, spline approximation methods and AV yields finite-order rates of 
convergence, which is observed in our calculations for S 3 . 

In order to approximate the initial data, or, in the case when the system 
involves distributed delays, the expansion coefficients of certain functions 
need to be computed. For AV and spline approximations such computations are 
relatively easy, since each element has local support. In our calculations 
for S 3 we used a Gauss quadrature rule [7]. In contrast, the Legendre 
polynomials are supported on the whole interval [- 2 , 0 ], and, for k large, 
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P k (x) is a rapidly oscillating function. A feasible algorithm for computing 
Legendre coefficients will be discussed in a forthcoming paper. However for 
AV and spline approximations, the expansion coefficients must be recomputed 
if N is changed. For the Legendre-tau approximation we can use the values 
which have already been computed for smaller N. 

N 

In the tables below we will use the following notation. 6 Is one of 
the differences 


| I (aJV(t) - x : (t)|, j=l, 2 , •••,!», 

k =0 J J 

N f-Vi N 

where a. (t) is the k ctl segment of dimension n in the solution vector a 

of the approximating system of ODEs (4.6), and x(t) is the true solution. 

Similarly, <5^ denotes one of the differences 
S 3 

lO N (0)wp (t) -Xj(t)|, j = l » 2, • • • ,n, 

where w N (t) is the coordinate vector of the approximate solution obtained by 
the cubic spline approximation method (S 3 ) (see, for the definition of 
S N (0), [2], p. 507). Note that dim w^ = n(N+3) and dim = n(N+l). 

Example 1 (Banks-Kappel [2], Example 1) 

In this example we study the equation for a damped oscillator with 
delayed restoring force and constant external force, 

d 2 d 

— 2 x (t) + dt + “ 10 » 


with initial conditions 
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x(8) = cos 9, — x(0) 

at 


- sin 9 for 0 e [-1,0]. 


Rewriting the above equation as a first order system we have 


d_ 

dt 


Xj (t) 


1 

o 

►— 

J 


"x^tr 

x 2 (t) 


H 

1 

o 


x 2 (t) 

_ _ 


L _J 


— — 


x^t-1) 


x 2 (t-l) 



r- n 
0 

+ 



10 


_ _ 


where x^ (t) = x(t) and x 2 (t) = ^ (t). 

Table II and Table III show the numerical results for x(t) and 

— x(t). For this example SP achieves the same accuracy with smaller N 
dt 

that Sg does. For SP, it appears that the rate of convergence is infinite- 
order. 


Table II 


t 

x(t) 

6 s 

b 3 

<1 

b 3 

,16 

% 

SP 

,8 

6 SP 

.25 

1.2704759 

0.00007 

0.00022 

0.000024 

0.00011 

0.0000047 

.5 

1.9936737 

0.00171 

0.00021 

0.000022 

0.00028 

0.0000012 

.75 

3.0614837 

0.00247 

0.00019 

0.000017 

0.00011 

0.0000176 

1.0 

4.3927203 

0.00083 

0.00024 

0.000012 

0.00027 

0.0000106 

1.25 

5.9259310 

0.00082 

0.00032 

0.000016 

0.00020 

0.0000016 

1.5 

7.6000709 

0.00090 

0.00005 

0.000010 

0.00007 

0.0000009 

1.75 

9.3440157 

0.00054 

0.00025 

0.000005 

0.00006 

0.0000002 

2.0 

11.0833011 

0.00013 

0.00003 

0.000015 

0.00003 

0 . 
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Table III 


t 

x(t) 

5 s 

5 s 

*16 

% 

SP 

CM 
CO CO 
«o 

.25 

2.06969 

0.00946 

0.00205 

0.00048 

0.00001 

0.00001 

.5 

3.64428 

0.00792 

0.00021 

0.00186 

0.00171 

0.00025 

.75 

4.84445 

0.00073 

0.00200 

0.00045 

0.00335 

0.00014 

1.0 

5.76581 

0.00790 

0.00192 

0.00046 

0.00348 

0.00031 

1.25 

6.45956 

0.00625 

0.00119 

0.00039 

0.00029 

0.00008 

1.5 

6.45956 

0.00185 

0.00245 

0.00062 

0.00078 

0.00002 

1.75 

7.0159957 

0.00061 

0.00003 

0.00067 

0.00075 

0 . 

2.0 

6.8497211 

0.00180 

0.00200 

0.00059 

0.00019 

0 . 

Example 2 
Next 

([2], Example 4) 
we use an example due to 

Popov for 

a degenerate 

system where we have 

(1,-2, -1)' 

r x(t) = 0 

for t > 2 

and all 

initial data (n.<|>) e 

Z. The 


equation is 


d_ 

dt 



"o 

2 

— % 

0 


0 

0 

o' 

x(t) = 

0 

0 

-1 

x(t) + 

1 

0 

0 


0 

0 

0 


0 

2 

0 




— 


- 


- 


x(t-l). 


We choose the initial data: 


n = cold, 1,1) and <K9) = 0 e K 3 


for -1 < 6 < 0 
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Note that (n,<J>) t 0(A). 

In Tables IV, V, and VI we give the numerical results. Here the initial 
function is not in the subspace for either method. With respect to 

* 2 (t) and xg(t) we can see that S 3 gives slightly better approximations 
than SP for t < 1. This is because X 2 and X 3 have a jump discontinuity 
in the derivative at t = 1. For SP, the convergence of the approximate 
solutions to the asymptotic solution ( 2 , 0 , 2 ) is quite rapid and seems to be 
infinite-order. 


Table IV 


t 

xj (t) 

S 3 

6 s 

b 3 

<•16 

% 

SP 

6 8 

SP 

*16 

S SP 

0.2 

1.36 

0.01211 

0.00493 

0.00247 

0.00140 

0.00041 

0.00001 

0.4 

1.64 

0.00889 

0.00528 

0.00001 

0.00146 

0.00158 

0.00001 

0.6 

1.84 

0.00445 

0.00334 

0.00190 

0.00547 

0.00098 

0.00016 

0.8 

1.96 

0.00980 

0.00074 

0.00146 

0.00233 

0.00111 

0.00019 

1.0 

2.0 

0.00020 

0.00221 

0.00054 

0.00339 

0.00033 

0.00009 

1.2 

2.0 

0.00725 

0.00424 

0.00175 

0.00237 

0.00002 

0.00001 

1.4 

2.0 

0.00160 

0.00460 

0.00085 

0.00064 

0.00003 

3xl0 -7 

1.6 

2.0 

0.00527 

0.00406 

0.00130 

0.00036 

0.00003 

13xl0 -7 

1.8 

2.0 

0.00331 

0.00281 

0.00242 

0.00028 

0.00001 

5xl0 -7 

2.0 

2.0 

0.00256 

0.00117 

0.00130 

0.00004 

5xl0 " 6 

lxl 0“ 7 

2.2 

2.0 

0.00306 

0.00049 

0.00096 

0.00004 

IxlO " 6 

9xl0 " 9 

2.4 

2.0 

0.00531 

0.00178 

0.00214 

0.00003 

2 xl 0“ 8 

-9 

2 x 10 * 

2.6 

2.0 

0.00212 

0.00246 

0.00115 

5xl0 -6 

9xl0 " 8 

IxlO " 9 

2.8 

2.0 

0.00055 

0.00251 

0.00092 

3xl0 " 6 

5xl0 ~ 8 

3xl0 -9 

3.0 

2.0 

0.00117 

0.00200 

0.00192 

2 xl 0~ 6 

2 xl 0“ 8 

0. 
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Table V 


t 

X2(t:) 

j 

j 

s s 6 

3 

S 4 

SP 

6 8 

SP 

,16 

6 SP 

0.2 

0.8 

0.00947 

0.00536 

0.00230 

0.02139 

0.00517 

0.00295 

0.4 

0.6 

0.01206 

0.00451 

0.00034 

0.02201 

0.00650 

0.00302 

0.6 

0.4 

0.00526 

0.00260 

0.00202 

0.00021 

0.00126 

0.00119 

0.8 

0.2 

0.00153 

0.00015 

0.00139 

0.02790 

0.00240 

0.00039 

1.0 

0. 

0.02133 

0.01895 

0.00949 

0.04081 

0.02072 

0.01034 

1.2 

0. 

0.01505 

0.00168 

0.00010 

0.00135 

0.00265 

0.00053 

1.4 

0. 

0.001369 

0.00068 

0.00094 

0.00440 

0.00085 

0.00013 

1.6 

0. 

0.00105 

0.00078 

0.00209 

0.00071 

0.00027 

0.00003 

1.8 

0. 

0.01149 

0.00105 

0.00193 

0.00067 

0.00006 

3 X 10~ 6 

2.0 

0. 

0.00444 

0.00123 

0.00020 

0.00042 

4 * 10 -6 

4 X 10~ 7 

2.2 

0. 

0.00636 

0.00108 

0.00173 

0.00004 

5 xl 0 -6 

3 X 10 -7 

2.4 

0. 

0.00521 

0.00066 

0.00202 

0.00007 

3 X 10 -6 

8 xl 0 -8 

2.6 

0. 

0.00205 

0.00019 

0.00073 

0.00003 

l x 10 -6 

l x 10~ 8 

2.8 

0. 

0.00407 

0.00021 

0.00084 

0.00005 

3 X 10 -7 

3 X 10 -10 

3.0 

0. 

0.00039 

0.00051 

0.00119 

0.00005 

3 X 10 -7 

6 X 10~ 10 
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Table VI 


t 

X 3 < t ) 

5 s 

b 3 

{ s 

b 3 

*16 

% 

6 4 

SP 

6 8 

SP 

.16 

6 SP 

0.2 

1.0 

0.01131 

0.00517 

0.00235 

0.04191 

0.01044 

0.00592 

0.4 

1.0 

0.01278 

0.00471 

0.00022 

0.04592 

0.01267 

0.00600 

0.6 

1.0 

0.00628 

0.00373 

0.00202 

0.00799 

0.00181 

0.00224 

0.8 

1.0 

0.00472 

0.00100 

0.00164 

0.04989 

0.00105 

0.00105 

1.0 

1.0 

0.04871 

0.03399 

0.01807 

0.08354 

0.04269 

0.02070 

1.2 

1.36 

0.02598 

0.00193 

0.00149 

0.00955 

0.00603 

0.00091 

1.4 

1.64 

0.02266 

0.00337 

0.00228 

0.00595 

0.00250 

0.00014 

1.6 

1.84 

0.00015 

0.00167 

0.00248 

0.00443 

0.00143 

0.00007 

1.8 

1.96 

0.01932 

0.00086 

0.00153 

0.00314 

0.00117 

0.00016 

2.0 

2.0 

0.01007 

0.00266 

0.00049 

0.00295 

0.00057 

0.00008 

2.2 

2.0 

0.00901 

0.00349 

0.00223 

0.00328 

0.00002 

0.00001 

2.4 

2.0 

0.00127 

0.00375 

0.00198 

0.00098 

0.00007 

1 * 10" 6 

2.6 

2.0 

0.00071 

0.00321 

0.00069 

0.00052 

0.00005 

3 xl 0~ 6 

2.8 

2.0 

0.00714 

0.00211 

0.00049 

0.00050 

0.00002 

9 xl 0" 7 

3.0 

2.0 

0.00465 

0.00076 

0.00060 

0.00011 

0.00001 

lxlO -7 


Examp le 3 ([2], Example 5) 

We consider the scalar equation 


x(t) = 5x( t) + x(t-l). 


with initial function 
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x(0) =5 for 0 e [-1,0]. 


The numerical results for this example can be found in Table VII and 
VIII. Again, we observe the quickness of the convergence of the approximate 
SP solution in this example. 


Table VII 


t 

x(t) 

*i 

S 3 

b 3 

.16 

% 

•i 2 

S 3 

6 64 

S 3 

0.2 

15.309691 

0.03 

0.00486 

0.00028 

0.000621 

0.000066 

0.4 

43.334337 

0.05 

0.00403 

0.00229 

0.000243 

0.000077 

0.6 

119.513222 

0.26 

0.01267 

0.00080 

0.000056 

0.000107 

0.8 

326.588900 

0.78 

0.01900 

0.00136 

0.000357 

0.000079 

1.0 

889.478955 

2.39 

0.09568 

0.00206 

0.000012 

0.000123 

1.2 

2420.772761 

7.28 

0.25894 

0.00031 

0.000845 

0.000149 

1.4 

6588.865818 

22.02 

0.77234 

0.00316 

0.001704 

0.000036 

1.6 

17934.153211 

65.58 

0.71924 

0.01228 

0.004029 

0.000051 

1.8 

48815.256906 

194.12 

6.70537 

0.03713 

0.013520 

0.001419 

2.0 

132871.377933 

570.95 

19.58949 

0.10230 

0.039730 

0.004005 





39 


Table VIII 


t 

x(t) 

6 4 

SP 

,8 

6 SP 

file 

SP 

6 32 

SP 

,64 

6 SP 

0.2 

15.309691 

0.04 

0.00823 

0.00021 

0.000291 

0.000046 

0.4 

43.334337 

6*10~ 5 

0.00316 

0.00076 

0.000368 

0.000058 

0.6 

119.513222 

0.22 

0.01232 

0.00196 

0.000193 

0.000042 

0.8 

326.588900 

0.60 

0.01626 

0.00285 

0.000452 

0.000066 

1.0 

889.478955 

1.67 

0.01392 

0.00220 

0.000334 

0.000049 

1.2 

2420.772761 

5.22 

0.00216 

0.00007 

0.000032 

0.000003 

1.4 

6588.865818 

15.56 

0.00108 

0.00004 

2xl0“ 7 

4xl0“ 7 

1.6 

17934.153211 

46.07 

0.00075 

0.00003 

0.000002 

00 

1 

0 

tH 

X 

1.8 

48815.256906 

135.58 

0.00002 

0.00002 

0.000003 

2xl0 -7 

2.0 

132871.377933 

396.72 

0.00069 

0.00009 

0.000014 

0.000002 


Example 4 


Here we deal with the equation which has multiple point delays 


x(t) = x(t) + 2x(t - V 2 ) + x ( t — 1 ) 


with initial function 


x(9) = 1 for 0 e [-1,0]. 


The numerical results in Table IX show that both methods work equally 


well 
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Table IX 


t 

x(t) 

j 

4 

3 

*16 

% 

*4 

6 -SP 

,8 

6 SP 

*16 

6 SP 

0.2 

1.885611 

0.002409 

0.000311 

0.000086 

0.006696 

0.000193 

0.000183 

0.4 

2.967299 

0.001213 

0.000022 

0.000403 

0.024093 

0.000916 

0.000748 

0.6 

4.331245 

0.004842 

0.001300 

0.000565 

0.021832 

0.000727 

0.000472 

0.8 

6.342954 

0.010185 

0.004028 

0.000170 

0.006117 

0.001556 

0.000436 

1.0 

9.278242 

0.005778 

0.003518 

0.000356 

0.011526 

0.001310 

0.000254 

1.2 

13.563777 

0.011080 

0.001466 

0.000741 

0.008807 

0.000168 

0.000027 

1.4 

19.903791 

0.010571 

0.000014 

0.000361 

0.007846 

0.000388 

0.000024 

1.6 

29.212354 

0.002034 

0.000492 

0.000141 

0.004891 

0.000145 

0.000025 

1.8 

42.845032 

0.015807 

0.001791 

0.000343 

0-015039 

0.000077 

0.000002 

2.0 

62.841170 

0.001850 

0.001460 

0.000224 

0.021896 

0.000072 

0.000013 


Table X 


t 

x(t) 

*! 

S 3 

*16 

°S 

b 3 

*32 

% 

A 8 

6 SP 

*16 

6 SP 

* 32 

6 SP 

0.2 

0.304159 

0.000447 

0.000099 

0.000051 

0.001490 

0.000040 

0.000105 

0.4 

- 0.463775 

0.000060 

0.000124 

0.000024 

0.001717 

0.000151 

0.000089 

0.6 

- 1.173353 

0.000424 

0.000006 

0-000014 

0.003747 

0.000519 

0.000058 

0.8 

- 1.796286 

0.000113 

0.000018 

0.000043 

0.004891 

0.000624 

0.000178 

1.0 

- 2.307740 

0.001049 

0.000400 

0.000101 

0.003346 

0.000082 

0.000115 

1.2 

- 2.497698 

0.002418 

0.000151 

0.000138 

0.001002 

0.000014 

0.000048 

1.4 

- 2.121921 

0.001361 

0.000205 

0-000073 

0.000652 

0.000013 

0.000003 

1.6 

- 1.236591 

0.000108 

0.000266 

0.000022 

0.000533 

0.000027 

0.000002 

1.8 

0.051303 

0.001261 

0.000063 

0.000019 

0.000482 

0.000037 

0.000003 

2.0 

1.606226 

0.001956 

0.000227 

0.000014 

0.000693 

0.000092 

0.000002 
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Example 5 

This example Is used in [1] for the identification problem which has a 

distributed delay and a discontinuous forcing function. This equation is 

, 0 

— x(t) = -3x(t) - / x(t+0)d9 + u^t), 

with initial function 

x( 0 ) =1 for 0 e [- 1 , 0 ], 


where u # 1 (t) is defined by 


u .i (t) 


1 

0 


on [ 0 , . 1 ] 
otherwise 


In Table X we give the numerical results. The approximations by S 3 
seem to be slightly better initially than by SP. This is because the solution 
x t (0) has a discontinuous derivative with respect to 0 for t < 1.1. The 
superiority of SP on the interval where the solution is infinitely 
differentiable is again observed for t > 1 . 2 . 

From the discussion and numerical results presented here we can conclude 
the following. For the Legendre-tau approximation, the effort in Implementing 
the algorithm is as much as for AV. But the accuracy of approximation is as 
good as that which can be obtained using high order spline approximation. 
When the solution to a problem is infinitely differentiable, the rate of 
convergence is faster than any finite power of (1/N). In addition, from our 
calculations, our method appears to be about four times as fast as the cubic 
spline approximation method. These properties are much more evident for 
examples (not presented in this paper) where optimal control and the 
approximation of eigenvalues were considered. 
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